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ABSTRACT 



Until now, axisymmetric, a-disc models have been adopted for calculations of the 
chemical composition of protoplanetary discs. While this approach is reasonable for 
many discs, it is not appropriate when self-gravity is important. In this case, spiral 
waves and shocks cause temperature and density variations that affect the chemistry. 
We have adopted a dynamical model of a solar-mass star surrounded by a massive 
(0.39 M0), self-gravitating disc, similar to those that may be found around Class 
and early Class I protostars, in a study of disc chemistry. We find that for each of 
a number of species, e.g. H2O, adsorption and desorption dominate the changes in 
the gas-phase fractional abundance; because the desorption rates are very sensitive 
to temperature, maps of the emissions from such species should reveal the locations 
of shocks of varying strengths. The gas-phase fractional abundances of some other 
species, e.g. CS, are also affected by gas-phase reactions, particularly in warm shocked 
regions. We conclude that the dynamics of massive discs have a strong impact on how 
they appear when imaged in the emission lines of various molecular species. 

Key words: stars: pre-main-sequence, stars: circumstellar matter, protoplanetary 
discs, astrochemistry 



. 1 INTRODUCTION 



In t he last several year s 

ter (|Dutrev et al.ll2007l: iHennin 



the Plateau de Bu re Interferome- 
et al. | l2010l ) and the Sub- 



millimeter Array (|Qi et al.l bOOSl l have provided images of 



protoplanetary discs in an increasing variety of molecular 
line emissions. The field will be advanced further by the 
com pletion of the Subm illimeter Array Disc Imaging Sur- 
vey (|Oberg et al.l 12010 1 . It will be revolutionised through 
the use of the Atacama Large Millimetre/submillimetre Ar- 
ray (ALMA), which will allow much weaker lines to be in- 
vestigated, and thus will enabl e the ch aracteri sation of a 
more diverse ran^e o f processes (jGuilloteau fc Dutrevll2008l : 
ISemenovet"aIll2008l ). 

Ma ny authors ( e.g., lAikawa et al 

2OO2I: (Markwick et all 



Gorti fc Hollenbach 
Nomura k Millar! 



.^^; 

I200J: 



iBergin et al.l 
lllgner et al 



20051: IWillacv et alJ 120061: IWillacv 



ail 



Nomura et al.ll2009l: IWoitke et al.ll2009l: 



iFogel et al.1 I2OIII : iHeinzeller et al.1 120111 



1999 . 



2003 



2004 



sented computational results for the chemical composition 
of protoplanetary discs. Because a-disc models were 
adopted, these studies are relevant to protoplanetary 
discs that transport mass locally, both instantaneously 
and averaged in time, and by extension, to discs that 
possess global axisymmetry. For this reason, a-disc models 
may be inappropriate during early phases of gas accre- 
tion on to protoplanetary discs if discs become massive 
enough to trigger gas phase gravitational instabilities 
(GIs) (|Vorobvov fc Basul l2005l . I2OO6I . l2009l : iBolevI l2009l ). 
GIs produce dynamic nonaxisymmetric str ucture in the 
form of spiral waves leading to shocks (see IPurisen et al.l 
[20071, for a review) and possibly, under some condi tions, to 
bound protoplanetary fragment s (Bos s 2001: Mayer et al.1 
2004 I2OO7I: IStarnatellos et al.1 I2OO7I: iBolev et al.1 I2OIOI : 



Vorobvov fc Basul I2OIOI : iBolev fc Durisenll2010l ) Frai 



Walsh et al. 



have pre- 



t at ion condition s are still under debate (see, e.g.,|BosL 
Cai et al.1 ^2010': 'P odolak et all [2010'; M eru &: Bat^ 



iBossi k 



2007; 



2011 



Lodato fc Clarke 20111 1 but, for non- fragmenting cases, it is 



generally agreed that gravitational torques caused by spiral 
waves lead to significant mass transport (jGammid I2OOII : 
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Lodato fc Ric3 1200^ fMeiia et alJ l2005l: JBolev et alJ l2006l . 






2007al : ICossins et al. 1 120091 : iMichael et al.ll201lh . 

Such mass transport may have imphcations for disc 
driven outbursts and recent, successful models for FU Orio- 
nis outbursts have been developed for the assumption that 
GIs dominate mass transport outs ide disc radii of a few AU 
in_non- fragmenting discs (Armita ge et alJ I2OOII : IZhu et alJ 
[2009. 2010). Simulations (Zhu et aL2oTo) including layered 
accretion in the Dead Zone (Gammie 1996 ) as well as trans- 
port by GIs show that, when discs accrete from rapidly ro- 
tating cloud cores, disc masses can become comparable to 
those of their central stars. Simple a-disc models for disc 
evolution overlook the episodic heating induced by GI spi- 
ral shocks in massive discs (^Bolev fc DurisenI bOOSl V These 
shocks can cause desorption of volatiles from dust grains and 
can trigger gas-phase chemical reactions that would not oth- 
erwise occur. Both effects can produce observable chemical 
signatures of disc dynamics. 

In this paper, we present results from a chemical model 
of a massive, young protoplanetary disc in which GIs cause 
the formation of spiral waves. Section [2] outlines the physi- 
cal and chemical models we have utilised, and the assumed 
initial conditions. Section [3] contains results for the time- 
dependent fractional abundances within an individual fluid 
element and column density maps of the entire disc for dif- 
ferent species. A comparison of results from this model with 
those of other models of disc chemistry are also given. Fi- 
nally, Section 2] presents conclusions based on the results 
and comments on further avenues of research. 



2 METHODS 

2.1 The dynamical model 

We use a hydrodynamic simulation of a massive (0.39 Mq) 
protoplanetary disc as the basis for the physical input to 
our chemical model. Most of the mass in the disc initially 
extends from a radial distance r '-^ 7 to 50 au from the cen- 
tral, solar-mass protostar. The system represents an early 
Class I object, which would likely evolve into an F-typ e main 
seque nce star. The disc is modelled with chymera (JBolevi 
[20071), which solves the equations of hydrodynamics with 
self- gravity on a fixed, cylindrical Eulerian grid. The compu- 
tational domain has 256, 128 and 64 zones in r, and z, re- 
spectively, with a physical resolution of Ar = Az = 0.25 au. 
Mirror symmetry about the midplane is assumed. All grid 
boundaries have outflow conditions allowing mass to leave 
the global disc simulation, where the inner boundary is set 
to approximately 4 au, but no mass is allowed to enter. The 
central star is allowed to move freely, and t he equ a tion o f 
motion for it is integrated as d escribed in iBolevi (|2009[ ). 
We use the iBolev et al.l (|2007bl l equation of state with a 
fixed ortho-to-para ratio of hydrogen molecules of 3:1 (which 
is only relevant for the evolution of the physical model). 
The fractional mass abundances of hydrogen, helium and 
more massive elements are set to X — 0.73, Y — 0.25, and 
Z — 0.02, respectively. Because the hydrogen is almost en- 
tirely in H2, the mean molecular weight is 2.33 amu. For cool- 
ing, w e use the radiative cooling approximation described in 
iBolevI |2009). An incident radiation field on the disc is as- 
sumed to have a black body spectrum with a temperature 
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Figure 1. Column density of nuclei (see Eqn.[T]for the relation- 
ship between number density and mass density) in the disc viewed 
from above at the end of the simulation, t = 388 yr. 



varying as Tirr = 140(r/au)~°'^ + 10K to account for heating 
by the star. 

The initial disc model is prepared using the method de- 
scribed in Boley & Durisen (2008). An analytic disc model 
is first created, with a iToomrd (|l964r ) Q ^-^ 1 for most 
of the disc. The spi ral instability sets in around Q ^ 1.7 
([Durisen et al.l 120071 ). so the initial model represents condi- 
tions before a strong burst of gravitational instabilities, as 
might be expected during the earliest stages of disc evolu- 
tion. A flat Q profile with the irradiation law noted above 
requires the surface density to follow, roughly, S oc r~^"^^. 
This initial configuration is assumed to undergo Keplerian 
rotation and have an adiabatic index of 5/3. These assump- 
tions about the initial rotation and value of the adiabatic 
index can lead to large radial oscillations at the start of the 
simulation due to the associated deviations from equilib- 
rium. To avoid some of this behavior, we evolve the analytic 
model in chymera at low azimuthal resolution (8 zones) for 
500 yr, which is about two orbital periods near r ^ 45 au. 
When the disc is loaded on to the higher resolution grid, a 
5% cell-to-cell random noise is added to the density distribu- 
tion, which seeds the growth of nonaxisymmetric structure. 
This is considered to be the time t = in the initial disc. The 
disc is then run at higher azimuthal resolution for an addi- 
tional 388 yr. Figure [1] shows the column density of nuclei 
for the disc viewed from above at the end of the simulation. 

At the start of the fuh simulation (t = 0), 1000 fluid 
elements are randomly distributed throughout the disc, 
weighted by mass. These fluid elements are evolved along 
with the main simulatio n, and their thermal h istories are 
recorded, as described in lBolev fc DurisenI (2008). The divi- 
sion of the disc into more fluid element would have allowed 
the capture of rare events such as shocks that are atypically 
strong for the radii at which they occur. This would lead to 
slightly greater average gas phase abundances of species that 
are hard to desorb. However, this would not affect our overall 
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Figure 2. The final location of the fiuid elements used to sample 
the disc, at t = 388 yr. 



conclusions. By the end of the simulation, 5 fluid elements 
were lost through grid boundaries, leaving 995 elements with 
complete thermal histories. 963 of these were used to study 
the chemical evolution in the disc. Results for 16 of the 32 
remaining fluid elements were not used to generate chem- 
ical results because the jump, for each of those 16, in the 
physical conditions from those at the end of the dynamical 
calculation back to those at the start of the dynamical cal- 
culation created difficulties for the integration of the chem- 
ical rate equations. Results for the other 16 were not used 
because spiked extrema in them caused problems for the 
chemical integration. These extrema occurred at tempera- 
tures of around 30 K - 40 K. The 32 elements were randomly 
distributed at the end of the simulation, and had not passed 
through shocks near the end of the dynamical calculation, 
so their influence on the chemical results was negligible. Fig- 
ure [2] shows the distribution of the final locations of the fluid 
elements used in the chemical modelling. 

Figure [3] shows the maximum line-of-sight tem- 
perature and the mass-averaged temperature Tn = 
(JnTdz]/(Jndzj along the line-of-sight for a face-on 
view of the final disc. Here n is the number density of nuclei 
(see Eqn.[T)), T is the temperature and dz is an infinitesimal 
path length along the line-of-sight. The maximum and mass- 
weighted line-of-sight temperatures were calculated from the 
full hydrodynamic model results, rather than from the prop- 
erties of the massive fluid elements used in the chemical 
calculations, because the full hydrodynamic model results 
provide higher spatial resolution. 

The spiral structure is clearly seen in the maximum 
line-of-sight temperature map. The structure is also seen in 
the mass- averaged temperature map, but is somewhat less 
sharply defined. The difference in the maps is due to the 
high temperature region being somewhat, but not highly, 
limited in extent. It lies near the midplane and contains 
only a fraction of the disc mass. 
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Figure 3. Maximum line of sight temperature (top) and mass 
averaged temperature (bottom) of the disc at t = 388 yr. 



Figure m shows results, taken from the complete hydro- 
dynamical data, for the temperature and number density of 
nuclei in planes containing the rotation axis and either the 
North-South axis (y) or the East- West axis (x) from Fig. (2] 
The temperature is less than 150 K outside the inner 20 au 
of the disc, while within it, the temperature can reach up to 
350 K. The number density plot shows that the highest den- 
sity gas lies in a ring-like structure at approximately 10 au 
from the centre. 

The temperature and density histories of the fluid ele- 
ments shown in Fig. [2] form the physical input to the chem- 
ical model. Figure \5\ shows the evolution of the tempera- 
ture and density of one of the fluid elements. This element 
begins at 0.8 au above the midplane, approximately 30 au 
from the centre of the disc. It follows a nearly circular path 
that slowly increases in radius. After one- and- a- half orbits, 
or 270 yr, it encounters the first shock. Both shocks do not 
affect the circular motion, but instead cause the element to 
rise to approximately 1.8 au above the midplane. The results 
for this track provide a useful example, because the fluid el- 
ement experiences a relatively quiescent period for 270 yr, 
before encountering the shocks. These shocks rapidly raise 
the temperature to a maximum of 140 K and the number 
density to 10^"^ cm""^. The temperature and density appear 
to peak at the same time and drop together. This is due 
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Figure 4. Slices of tlie disc interior sliowing number densities and temperatures at t = 388 yr. Tlie x-axis liere corresponds to a slice 
along the a::-axis of Fig. [T] where y = 0, and vice versa. 
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Figure 5. Temperature and number density history of a fluid 
element from the disc. This particular element encounters a shock 
at about 270 and 350 yr. 



to the increase in the pressure caused by a shock driving 
an expansion of gas parallel to the shock and obliquely to 
the mid-plane of the disc that begins almos t immediately 
behind the shock fsee lBolev fc Durisenll2008l ). 



2.2 The chemical model 

The chemical computation is based on a network of rate 
equations involving 1334 reactions and 125 species contain- 
ing the elements H, He, C, O, N, Na and S. These reac- 
tions were originally selected from a subset of the UMIST 



Rate 95 database ([Millar et al.lll997l V Data from the Kinetic 
Database For Astro chemistry, KIDA[j, were used to update 
some of the rates and rate coefficients. Thermal desorption 
processes (discussed in detail below) and some three-body 
reactions (see Appendix [BJ were added to the network. 

At each time-step in the chemical network, the tem- 
perature and mass density of each fluid element are derived 
with cubic spline fits to the data provided for the element 
from the physical model. These interpolated values of the 
temperature and mass density are then used in the time- 
integration of the rate equations. The number density of 
nuclei, n, is calculated from the mass density, p, by 

n = — p = riH + nue + nz, (1) 

where nn, nue and nz are the number densities of hydro- 
gen, helium and more massive nuclei, respectively. We have 
assumed a molar mass of M = 1.28 g mol~^, which is ap- 
propriate for a relative abundance of helium-to-hydrogen of 
0.09 with trace amounts of more massive elem ents. Integra- 
tions were performed with the dvode package ([Brown et al.l 
Il989l ) to obtain the fractional abundance, X{i) — n{i)/n, of 
each species. 

The rate equation for the gas-phase fractional abun- 
dance of the ith species is 



d_ 
dt 



X{2) 



k{j)X(l)X{m)n - Y^ k{j')X{i)X{m) 

j' ,m 

2j2k{f)X{ifn - T,r{i)X{{) + Ssii) + S'a,d(i 






(2) 



where k{j) is the rate coefficient of the jth reaction and the 
summations are restricted so that only reactions involved 

^ jhttp : //kida . obs . u-bordeauxl . f r | 
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in the formation or removal of the ith species are included. 
rcr(0 is the rate at which direct cosmic-ray ionisation and 
cosmic-ray-induced photoemission remove species i. Ssii) 
and 5'a,d(0 ^^^ source terms due to three-body reactions and 
adsorption on to and desorption from grain surfaces, respec- 
tively. We assume that the only third body of importance 
is H2, an assumption which introduces only an insignificant 
error. Thus, 

Ssii) - Yl HJ)X{l)X{m)X{R2)n (3) 

j,l,m 

- Y^ k(f)X{i)X{m)X{ll2)n 

j' ,m 

2.2.1 Rate coefficients and rates 

Here we present a brief summary of the forms of the gas- 
phase rate coefficients and rates that we have adopted. More 
detail can be foun d in the primary p aper on the UMIST 
Rate 95 database (JMillar et al.l 119971 ) . The reaction rates 
and coefficients that we use are provided onlincl- 

The rate coefficient for the jth two- or three-body re- 
action is of the standard Arrhenius form 



kU) 



<J) ( 



— ) 
3007 



^U) 



exp 



-7(i) 



(4) 



where a{j) is the room temperature rate coefficient of the 
reaction (at 300 K), I3(j) describes the temperature depen- 
dence and 7(j) is the activation energy of the reaction. 

The rate for the destruction of species i due to cosmic 
rays is given by 

CP(i) 



TcrW =a(i)C + 



1-A' 



(5) 



where a(i) is a proportionality constant, ( = 10~^^s~^ is 
the cosmic ray ionisation rate, P(i) is a constant and A is 
the dust grain albedo, which we take to be 0.5. Because 
the two terms on the right represent the physically distinct 
processes of direct cosmic ray induced ionisation and of de- 
struction due to photoemission induced by the collisions of 
molecular hydrogen with energetic electrons produced by 
the ionisation, data tables give a{i) and P{i) separately. 

Though photoabsorption of radiation from external 
sources will affect the chemistry in the outer layers of a disc, 
we focus on the bulk of the disc material and assume that 
it is well shielded from external sources of photons. Thus, 
we neglect photoabsorption, other than that of cosmic-ray- 
induced photoemission. 

2.2.2 Gas-grain reactions 

With minor exceptions, we assume that no surface chem- 
istry occurs, though we recognise that surface chemistry 
is important in establishing the chemical initial condi- 
tions that we adopt. The exceptions concern the neutrali- 
sation of ions, which we assume occurs when they are ad- 
sorbed. While studies of the chemistry of discs including 



surface reactions and grain evolution have been performed 
(ISchrever et al. 20081: iHenning et aL 2010l: ISemenov et al.l 



'2010': 'Vas vunin et aljEoill : ISemenov fc Wiebell201lh . such 
chemistry remains very uncertain. This is why we have as- 
sumed that the initial abundances reflect those of cometary 
ices (see Section [2. 3p . 

Sa,d{i) has contributions due to adsorption and desorp- 
tion, and we take 



^a,d(i) = ^d(i)-&(i), 

where 



Sa — 7va'^S(i)r]\ X(i)n. 



(6) 



(7) 



The sticking coefficient S{i) is set to unity, and we take the 
ratio of the number density of dust grains to the number 
density of nuclei to be 77 = 3.3 x 10~^^, which is appropriate 
for a dust grain radius and grain-to-gas mass ratio of about 
a = 0.1 /im and 0.01, respectively. Tg is the gas temperature, 
k is the Boltzmann constant, /x is the molecular mass in amu 
and mu is the mass of a hydrogen atom. 

We tre a t the rmal desorption in the same way as 
IVisser et aJ] (|2009l ). and 



Sd{i) = 1.26 X 10"^^(cr/cm"^)/(i)zyo(i) exp 



-E^(i) 
kTd 



,(8) 



where a is surface density of binding sites, which we take 
to be 1.5 X lO^^cm"^. ^b(0 is the binding energy of the 
ith species on the surface of the dust grain, Td is the dust 
temperature (which we assume to be in equilibrium with 
Tg), the factor f{i) represents the fraction of the surface of 
the dust grain covered by the ith species, given by 



f{^) 






(9) 



where X^{i) is the solid fractional abundance of the ith 
species and A^b is the typical number of binding sites per 
grain, taken to be 10^. The characteristic vibrational fre- 
quency of the species attached to the grain is given by 



i^o[i) 



2aE^{i) 



(10) 



where m{i) is the mass of the ith species fHasega wa et al.l 
1992). The b in ding energies were taken from 
iHollenbach et al.l ((20091), references therein, and the 
OSU databas4!l. 



2.3 Chemical initial conditions 

The initial gas-phase fractional abundances were assumed 
to have the same ratios to one another as the fractional 
abundances of the corresponding; ices in co met Hale-Bopp 
as given by lEhrenfreund fc CharnlevI (|2000l l , and they are 
given in Table [TJ Because comets are thought to have un- 
dergone significant chemical processing at the edges of discs, 
this assumption may not be entirely appropriate. However, 
comparisons between the compositions of cometary ices and 
interstellar ices imply a general consistency between the two. 



' |http : //www . ast . leeds . ac . uk/~ pyj di/chemdisc /network . txt | 



|http : //www . physics . ohio- state . edu/~ eric/research . html | 
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Table 1. Initial fractional abundances {X{i) = n(i)/nii). Note 
that a(b) = a x 10^ . 



Species 


Abundance 


Species 


Abundance 


He 


l.OO(-l) 


H2CO 


1.83(-6) 


C 


3.75(-4) 


CO2 


3.67(-5) 


CO 


3.66(-5) 


HCN 


4.59(-7) 


CH4 


1.10(-6) 


HNC 


7.34(-8) 


N 


1.15(-4) 


S 


1.62(-5) 


NHs 


3.30(-6) 


H2S 


2.75(-6) 


O 


6.74(-4) 


SO 


1.47(-6) 


H2O 


1.83(-4) 


S02 


1.84(-7) 


Na 


3.50(-5) 


ocs 


3.30(-6) 




e discrepancies exi 






though som 


St (lEhrenfreund & CharnlevI 


|2000|;|Ehrenfreund & Schutte||200(]D . 





2.4 Timescales 

As mentioned above, the hydrodynamical simulation follow- 
ing the disc evolution is run at higher resolution for 388 yr. 
This is far longer than the orbital period of about 4 yr at the 
inner boundary and somewhat shorter than the orbital pe- 
riod of 390 yr at 60 au. Nearly all of the mass is at less than 
50 au where the period is about 300 yr, and the gravitational 
instability develops fully and leads to well established spi- 
ral structure in about 200 yr. The spiral structure roughly 
co-rotates with the material between about 30 and 40 au; 
the orbital period at 35 au is 180 yr. The disc is not steady, 
and the mass infall rate varies with radius and time from a 
factor of a few smaller than to a factor of a few larger than 
IO^^Mq, which implies a radial flow time of the order of 
10^ yr. These timescales are much longer than the shortest 
chemical timescales which are associated with adsorption 
and desorption; using Eqn. 7, one finds that the adsorp- 
tion timescale for a number density of 10^"^ cm""^ is roughly 
an hour. The desorption timescale is very much larger or 
smaller, depending on the temperature. 

388 yr is sufficient to allow the spiral shocks to develop 
fully. This disc is already exhibiting spiral structure that is 
typical in Gl-active discs over many orbits. The purpose of 
the hydrodynamical simulation is to provide self-consistent 
shock profiles for the chemical model, which was achieved in 
a relatively short period of evolution. 

We follow the chemistry of each fluid element through 
ten cycles of the temperature and number density history 
associated with it (e.g., the history shown in Fig. O, where 
in each subsequent cycle, we used the final fractional abun- 
dances from the previous cycle as initial input. Of course, 
this resulted in rapid variations in the physical conditions, 
as each new cycle began. However, we followed the chemi- 
cal evolution in this way to keep the amount of data pro- 
duced by the model at manageable levels. We found that 
the chemical behaviour for a fluid element became periodic 
by the time that the chemistry had gone through ten cycles. 
All results presented here are from the final cycle, and t = 
occurs at the beginning of that cycle. 

At the start of each integration, we used a logarith- 
mically increasing time-step, which was initially 90s. This 
allowed initially rapid reactions, such as the adsorption of 
species on to grains at the beginning of the first cycle, to be 



followed with sufficient resolution. We limited the maximum 
time-step to approximately 10^ s to avoid missing details in 
rapidly changing shock features. 



3 RESULTS 

The main results of the paper are column density maps for a 
variety of gas-phase species. To obtain them, we followed the 
chemical evolution of each fluid element. Before presenting 
the column density maps, we consider the chemistry in one 
fluid element as an illustrative example of the effect of shocks 
on the chemistry of the material. 



3.1 An individual fluid element 

Figure [6] shows the fractional abundances of 17 species (CO, 
SO, NH3, H2O, H2S, ocs, O2, HCO, HCO+, HNO, SO2, 
CS, HCS, HCS+, HCN, HNC and OCN) as functions of 
time during the final cycle for the fluid element for which 
the density and temperature are given as functions of time 
in Fig. \5\ The shocks induce increases in the gas-phase frac- 
tional abundances of all 17 species. However, the increase in 
X(CO) is too small to be apparent from the figure, because 
nearly all CO is in the gas-phase even in the coldest regions 
of the disc. 

The increases of the fractional abundances of species 
for which results are shown in the upper two right hand 
panels are due almost entirely to thermal desorption from 
the surfaces of dust grains. We define Ng{i) to be the av- 
erage number of molecules of the zth species on a grain. 
For each species for which results are given in those panels, 
Ng{i)r] + X{i) is nearly constant. This is due to adsorption 
and desorption being the only processes significantly affect- 
ing the gas-phase fractional abundances of these species, and 
the fractional abundance of each of them depends almost en- 
tirely on the temperature only. 

The fractional abundances of species for which results 
are shown in the other panels are affected significantly by 
gas-phase reactions, as well as desorption and adsorption. 
The higher temperatures in the shocks allow some reactions 
with activation energies and some endothermic reactions to 
occur at significant rates. Also species which undergo reac- 
tions that are exothermic and have small or no barriers are 
desorbed into the gas-phase in shocked regions. 

HCS is one species having a fractional abundance that is 
affected by gas phase reactions. It is formed by the reaction 
of S and CH2. CS is formed by reactions of C with SO, 
and S with CN in almost equal proportions. OCN is formed 
by the reaction of CN and O2. HNO forms by the reaction 
of NH2 and O, and HCO forms mainly via NH^ reacting 
with H2CO. The reaction of O with OH contributes to O2 
production throughout much of the disc; O2 is also formed 
in the shocked regions by the reactions of He^ with CO2 
and SO2. 



3.2 The entire disc 

A major aim of the work reported here is the identification of 
the general features that should appear in images of gravita- 
tionally unstable protoplanetary discs obtained in molecular 
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Figure 6. Fractional abundances of species as functions of time for the fluid element for which the physical structure is shown in Fig. 
O The relatively cold, quiescent period for the flrst 270 yr shows the abundance of most species decreasing due to adsorption, or staying 
constant. 



line emissions. Thus, we present model column density maps 
for a number of molecular species. 

The final spatial positions of the fluid elements consti- 
tute an irregular grid for which the model fractional abun- 
dances are known. The QHULL and QGRId3 interpolation rou- 



tines contained within the IDL 5.5 libraries were used to 
obtain the fractional abundances at the points of a three- 
dimensional regular Cartesian grid. The mass density at 
each of those grid points was calculated from the mass den- 
sity distribution at the final time given by the full hydrody- 
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namic results. We used the mass density distribution given 
by the full fluid results because the grid for which it was 
calculated is much more uniform than the grid for which 
the chemical results are available. Thus, for regions of low 
mass density, we were able to obtain more accurate inter- 
polated results for the mass density than for the fractional 
abundances. 

The column density N{i,x,y) of the ith species on a 
line of sight perpendicular to the disc is 



Table 2. Maximum abundances reported in this work, W2010 
and 12004. A dash signifies no data for that species were available 
in the original paper. 



N{i,x,y) 



I 



nX(i, x,y, z) dz. 



(11) 



We assumed that n and all X{i^x^y^z)s are even functions 
of z when performing the integrals. Figures [71 and [8] show 
the column densities of selected species. For each species for 
which results are given in Fig. [71 Ng(i)r] + X(i) is nearly 
constant. This is due to adsorption and desorption being 
the only processes significantly aflFecting the gas-phase frac- 
tional abundances of these species. Thus, the structure seen 
in each of these column densities is due to the variations in 
temperature and density and is not influenced by gas-phase 
chemistry. 

Ng{i)r]-\- X{i) is not constant for each of the species for 
which Fig. [8] shows a gas-phase column density map. Gas- 
phase chemistry significantly aflFects the gas-phase fractional 
abundance of each of these species. We have already given 
some explanation for the evolution of most of the species 
for which results are displayed in Fig. [8] in Section 13.11 The 
HCO^ is anticorrelated with the H2O and other species with 
which it charge exchanges. HCS^ is also removed by charge 
exchange with H2O, but its distribution diflFers somewhat 
from that of HCO^. This is due to CS being abundant in the 
gas-phase only in regions in which it is eflficiently desorbed 
and CO being abundant everywhere in the disc. Reactions 
of HJ with CS and CO produce these ions. 

3.3 The effect of three-body reactions 

The chemistry in each of several fluid elements was calcu- 
lated both with the three-body reactions included and with 
them excluded. As expected, the largest diflFerences occur in 
the denser inner regions of the disc and spiral arms. Even 
there, the eflFects were insignificant for all species except 
some with low fractional abundances {X(i) < 10"^*^). 

A reduction of the fractional abundance of O2 in one 
element from approximately 10~^^ to 10~^^ in a shocked 
region and an enhancement in the fractional abundance of 
NH from approximately 10~^^ to 10~^^ were amongst the 
most pronounced eflFects of including the three-body reac- 
tions. Species with fractional abundances higher than ap- 
proximately 10~^° were not noticeably aflFected by the addi- 
tion of the three-body reactions. The results that we present 
were obtained with the three-body reactions included. 

3.4 Comparison with other models 

The three-dimensional nature of our chemical model and 
the qualitatively diflFerent dynamics of the disc make a direct 
comparison of our results with the previously existing chem- 
ical results for axisymmetric models of lower mass discs diflfi- 
cult. However, we give a rough comparison between the max- 
imum fractional abundances of species given by our model 
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and some obtained by others who have investigated chem- 
istry for disc models. The results taken from the other pa- 
pers are somewhat inaccurate because they were obtained 
through the examination of flgures from which precise in- 
formation was hard to infer. Table [2] shows the maximum 
fractional abundances recorded at the end of our simula - 
tion and the co rresponding values from Walsh et al.l (|2010h 
and llkner e~ l. (2004) (hereafter W2010 and 12004 respec- 
tively). Our results are comparable to those of the other 
models for NH3, SO2, H2O, CO2, CN, CS and HCN. We 
obtain a higher fractional abundance of HCO^ than 12004, 
but less than reported in W2010. Our fractional abundances 
of C2H and N2H^ are lower than those of W2010. However, 
we report a much higher fractional abundance of H2CO. 

Some of the discrepancies between the results of W2010 
and 12004 arise because W2010 considered a much more ex- 
tensive fraction of the disc than 12004 who studied only the 
material at radii less than 10 au. We studied material at radii 
between 10 and 60 au, where spiral structure and shocks are 
important in our dynamical model. Thus, in further dis- 
cussion of the diflFerences between our results and those of 
others, we focus on a comparison of our results and those of 
W2010. 

Some of these diflFerences are due to the assumptions 
that we made concerning the initial chemical conditions. As 
mentioned earlier, in our model a number of species have 
gas-phase abundances that depend only on the assumed ini- 
tial conditions and the balance between desorption and ab- 
sorption. H2CO in particular is such a species. So its high 
abundance in our model compared to that found by W2010 
is due to our assumption of a higher H2CO abundance and 
neglect of species as massive as CH3OH. Similarly, our as- 
sumption that no nitrogen is initially in N2 leads to a lower 
abundance of N2H+ than W2010 found. 

The high values of the HCO+ and C2H given for W2010 
occur in regions with densities that are about four orders of 
magnitude lower than any that we consider; these regions 
have correspondingly lower total column densities than the 
much denser regions do and in the W2010 model, the chem- 
istry in them is aflFected by the diflFusion that W2010 assume 
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Figure 7. gas-phase column densities of molecules having gas-phase fractional abundances determined primarily by desorption and 
adsorption, at t = 388 yr. Distances from the disc centre are in au. 



to occur, but which we do not. The mid-plane value obtained 
by W2010 for the fractional abundance of HCO^ for radii 
of 10 to 60 au is comparable to the maximum value that we 
report for it; this is to be expected because in both models 
most of the CO is in the gas phase, and HCO^ is formed 
by a simple ionisation triggered reaction sequence. The mid- 
plane value obtained by W2010 for the fractional abundance 
of C2H for radii of 10 to 60 au is less than the maximum value 
that we report for it; this is due to the desorption, caused 
by shock induced heating, occurring near the mid-plane in 
our models. 



4 CONCLUSIONS & FUTURE WORK 

We have constructed a chemical model of a 0.39 M© proto- 
planetary disc, surrounding a solar-mass star, representative 
of a Class or early Class I system. In the disc, gravitational 
instabilities cause spiral waves. The shocks associated with 
these waves induce the desorption of various chemical species 
from the surface of dust grains and increase their gas-phase 
abundances. Though the gas-phase fractional abundances 



of some of the species are not significantly affected by gas- 
phase chemistry, some of the desorbed species are reactive. 
In some cases, the elevated temperatures in the shocked re- 
gions enhance the reactivity. 

Because most of the mass is concentrated in spiral struc- 
tures, all of the gas-phase molecular column densities show 
spiral structures. However, the structures are limited in ex- 
tent in some species, e. g. H2O, which possess the highest 
binding energies to grains, therefore higher temperatures are 
needed for desorption to occur. This also limits the extent of 
structure in species which form in the gas-phase from these 
tightly bound species. Consequently, maps of the emissions 
of a number of species will reveal where shocks of differ- 
ing strengths are and, thus, serve as diagnostics of the disc 
dynamics. 

We find that three-body reactions have little effect on 
the chemistry of species with fractional abundances above 
-j^Q-10 xhough they are most important in the hotter, dens- 
est regions of the disc, they do not alter the overall chemistry 
of the disc significantly. 

A direct comparison of our results with those of other 
models is not straightforward, due largely to the very differ- 
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Figure 8. Column densities in the disc for several gas-phase molecules whose total abundance varies with time, at t = 388 yr. Distances 
from the grid centre are in au. 



ent nature of the disc dynamics that we have used. Some dif- 
ferences between the peak gas phase fractional abundances 
that we obtained and those reported in other papers are due 
to the assumptions that we have made about chemical initial 
conditions. However, the major differences in the chemistry 
away from the outer boundaries of the disc arise from the 
role that some other modellers assume that diffusion plays 
in enhan cing the richness of the chemistry in the disc inte- 
rior fe.g. Illgner et al.ll2004l : iHeinzeller et al.ll201l! l: in many 
of these models, the chemistry at the midplane is not very 
rich over a very large fraction of the disc. In contrast, our 
model gives rise to substantial gas-phase abundances near 
the midplane, even though we have not assumed that ef- 
ficient microscopic mixing is a consequence of large-scale 
turbulence. Though inclusion of such mixing would be dif- 
ficult in an approach in which individual fluid elements are 
followed, we have not neglected that type of mixing for that 
reason. Rather, we have chosen to focus on how gravita- 
tional instability generates shocks which induce a rich gas- 



phase chemical composition in the disc interior, even near 
the midplane (see Appendix lAJ) . 

ALMA will revolutionise observational studies of discs. 
If used to observe a disc in the Taurus- Auriga cloud complex, 
ALMA, with its 5 milliarcsec resolution, will allow the map- 
ping of features on scales of about 1 au, which is smaller than 
the widths of the prominent spiral structures of our model. 
We plan to use our results in radiative transfer calculations 
to obtain reliable estimates of the detectability of structure 
in different species. 
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APPENDIX A: OTHER PERSPECTIVES 

Figure lAll shows the column densities of selected species in 
the disc as viewed along the ^/-axis, toward 2/ = 0, in Fig. [1] 
The majority of the edge-on maps show little structure (as 
with CO and H2O), though some, such as the H2CO and 
SO maps, show slightly higher column densities on the right 
side due to the presence of a large spiral arm. 

Figure IA2I shows slices of the fractional abundances in 
the interior of the disc for HON and H2CO. These show 
little structure in the disc interior when compared with 
other models of axisymmetric discs. These shock-desorbed 
molecules show little structure in the disc interior compared 
with what is seen in axisymetric a-disc models. 



APPENDIX B: THREE-BODY REACTIONS 

Table lBTI lists the three-body reactions included in t he chem- 
ical network and taken from I Willacy et al.l (|l998h and the 
more recent UMIST Database for Astrochemistry, UDfAJj. 
The subset includes only the reactions possessing activation 
energies of less than 10^ K. 



This paper has been typeset from a T^]X/ MJ]X file prepared 
by the author. 
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Figure Al. Column densities of the disc as viewed along the y-axis of Fig. [T] towards y = 0, at t = 388 
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Figure A2. Slices of the fractional abundance of disc interior, oriented as in Figure lU interpolated from the positions of the fluid 
elements at t = 388 yr. 
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Table Bl. Three-body reactions incl uded in the chemical networ k. In all but one of t he reactions the third body is assumed to be 
molecular hydrogen. References are W: IWillacv et al.l (|l998l ) and U: IWoodall et al.l (|2007l ). 
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